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Abstract 

Many problems in finance require the information on the first passage time (FPT) 
of a stochastic process. Mathematically, such problems are often reduced to the 
evaluation of the probability density of the time for such a process to cross a certain 
level, a boundary, or to enter a certain region. While in other areas of applications 
the FPT problem can often be solved analytically, in finance we usually have to 
resort to the application of numerical procedures, in particular when we deal with 
jump-diffusion stochastic processes (JDP). In this paper, we propose a Monte- 
Carlo-based methodology for the solution of the first passage time problem in the 
context of multivariate (and correlated) jump-diffusion processes. The developed 
technique provide an efficient tool for a number of applications, including credit 
risk and option pricing. We demonstrate its applicability to the analysis of the 
default rates and default correlations of several different, but correlated firms via a 
set of empirical data. 

Keywords: First passage time; Monte Carlo simulation; Multivariate jump-diffusion 
processes; Credit risk 



1 Introduction 



Credit risk can be defined as the possibility of a loss occurring due to the financial 
failure to meet contractual debt obligations. This is one of the measures of the like- 
lihood that a party will default o n a financial agreement. There exist two classes of 
credit risk models dZhoull2001bl) . structural models and reduced form models. Struc- 
tural models can be traced back to the influential works by Black, Scholes and Merton 
dBlacketalll20(51 iMertonl, |2005l). w hile reduced form models seem to originate from 
contributions by J arrow et all ( 1997 ). The major focus in this contribution is given to 
structural credit-risk models. 
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In structural credit-risk models, a default occurs when a company cannot meet its 
financial obligations, or in other words, when the firm's value falls below a certain 
threshold. One of the major problems in the credit risk analysis is when a default 
occurs within a given time horizon and what is the default rate during such a time 
horizon. This problem can be reduced to a first passage time (FPT) problem, that 
can be formulated mathematically as a certain stochastic differential equation (SDE). 
It concerns with the estimation of the probability density of the time for a random 
process to cross a specified threshold level. 

An important phenomenon that we account for in our discussion lies with the fact 
that, in the market economy, individual companies are inevitably linked together via 
dynamically changing economic conditions. Therefore, the default eve nts of compa- 
nies ar e often correlated, especially in the same industry. IZhoul (1200 1 aft and iHull et al 
(1200 ll) were the first to incorporate default correlation into the Black-Cox firs t pas- 
sage str u ctural model, bu t they have not included the jumps. As pointed out in IZhou 
d2001bh : IKou et al l d2003l) . the standard Brownian motion model for market behavior 
falls short of explaining empirical observations of market returns and their underlying 
derivative prices. In the meantime, jump-diffusion processes (JDPs) have established 
thems elves as a sound alternative to the standard Brownian motion model ( lAtiva et all 
20051) . Multivariate jump-diffusion models provide a convenient framework for in- 



vestigating default correlation with jumps and become more readily accepted in the 
financial world as an efficient modeling tool. 

However, as soon as jumps are incorporated in the model, except for very basic 
applications where analytical solutions are available, for most practical cases we have 
to resort to numerical procedures. Examples of known analytical solutions include 
problems where the jump sizes are doubly exponential or exponentially distributed 
(IKou et all |2003) as well as the jumps can have only nonneg ative values (assu ming 
that the crossing boundary is below the process starting value) dBlake et allll973l) . For 
other situations, Monte Carlo methods remain a primary candidate for applications. 

The conventional Monte Carlo methods are very straightforward to implement. We 
discretize the time horiz on into N intervals with N being large enough in order to 
avoid discretization bias dKloeden et all 120031) . The main drawback of this procedure 
is that we need to evaluate the processes at each discretized time which is very time- 
consuming. Many researchers have contributed to the field of enhancemen t of the ef- 
ficiency of Monte Carlo simulations. Among others, iKuchler et all (120021) discussed 
the solution of SD Es in the framework of weak discrete time approximations and 
Liberati et a 3 d2006l) considered the strong approxim ation where the SDE is driven by a 



high i ntensity Poisson process. Atiya and Metwally dAtiya et a l l2005llMetwallvetall 



2002) have recently developed a fast Monte Carlo-type numerical methods to solve the 



FPT problem. In our recent contribution, we reported an extension of this fast Monte- 
Carlo-type method in the context of multiple non-correlated jump-diffusion processes 
dZhanget alll2006l) . 

In this contribution, we generalize our previous fast Monte-Carlo method (for non- 
correlated jump-diffusion cases) to multivariate (and correlated) jump-diffusion pro- 
cesses. The developed technique provides an efficient tool for a number of applica- 
tions, including credit risk and option pricing. We demonstrate the applicability of this 
technique to the analysis of the default rates and default correlations of several different 
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correlated firms via a set of empirical data. 

The paper is organized as follows, Section [2] provides details of our model in the 
context of multivariate jump-diffusion processes. The algorithms and the calibration of 
the model are presented in Section[3] Section[4]demonstrates how the model works via 
analyzing the credit risk of multi-correlated firms. Conclusions are given in Section[5] 



2 Model description 

As mentioned in the introduction, when we deal with jump-diffusion stochastic pro- 
cesses, we usually have to resort to the application of numerical procedures. Although 
Monte Carlo procedure provide a natural in such case, the one-dimensional fast Monte- 
Carlo method canno t be directly gener alized to the multivariate and correlated jump- 
diffusion case (e.g. IZhang et all d2006l) ). The difficulties arise from the fact that the 
multiple processes as well as their first passage times are indeed correlated, so the 
simulation must reflect the correlations of first passage times. In this contribution, we 
propose a solution to circumvent these difficulties by combining the fast Monte-Carlo 
method of one-dimensional jump-diffusion processes and the generation of correlated 
multidimensional variates. This approach generalizes our previous results for the non- 
correlated jump-diffusion case to multivariate and correlated jump-diffusion processes. 

In this section, first, we present a probabilistic description of default events and 
default correlations. Next, we describe the multivariate jump-diffusion processes and 
provide details on the first passage time distribution under the one-dimensional Brow- 
nian bridge (the sum-of-uniforms method which is used to generate correlated mul- 
tidimensional variates will be described in Section [XTT i. Finally, we presents kernel 
estimation in the context of our problem that can be used to represent the first passage 
time density function. 



2.1 Default and default correlation 

In a structural model, a firm i defaults when it can not meet its financial obligations, or 
in other words, when the firm assets value Vi(t) falls below a threshold level Dy^t). 
Generally speaking, finding the threshold level Dy i (t) is one of the challenges in us- 
ing the structural methodology in the credit risk modeling, since in reality firms often 
rearrange their liability structure when they have credit problems. In this contribu- 
tion, we use a n exponential form defining the threshold level Dy^t) = Kiexp^t) 



as proposed by Zhou (2001a,), where 7^ can be interpreted as the growth rate of firm's 
liabilities. Coefficient Kj captures the liability structure of the firm and is usually de- 
fined as a firm's short-term liability plus 50% of the firm's long-term liability. If we 
set Xi(t) = ln[Vj(i)], then the threshold of Xi(t) is D. t (t) = + ln(Kj). Our main 
interest is in the process Xi(t). 

Prior to moving further, we define a default correlation that measures the strength of 
the default relationship between different firms. Take two firms i and j as an example, 
whose probabilities of default are P, and Pj, respectively. Then the default correlation 



3 



can be defined as 



p. 



Pij 



(1) 



Pi). Let us assume 



^(1-^)^(1-^) 
where Py is the probability of joint default. 

From Eq. © we have = PiPj + ,»,, x ~\ /',:/',! 1 
that Pi = Pj = 5%. If these two firms are independent, i.e., the default correlation 
Pij = 0, then the probability of joint default is P^ = 0.25%. If the two firms are 
positively correlated, for example, p^j = 0.4, then the probability of both firms default 
becomes P^ = 2.15% that is almost 10 times higher than in the former case. Thus, the 
default correlation py plays a k ey role in the j oint default with impo rtant implications 
in the field of credit analysis. IZhoul ( 12001 al) and Hull et all (120011) were the first to 
inc orporate defau lt correlation into the Black-Cox first passage structural model. 



Zhoul d2001al) has proposed a first passage time model to describe default correla- 



tions of two firms under the "bivariate diffusion process": 







Pi 


dt + n 




X 2 (t) 






dz 2 



(2) 



where /ii and [i 2 are constant drift terms, z\ and z 2 are two independent standard 
Brownian motions, and O is a constant 2x2 matrix such that 



n ■ fi' = 



p(J\(J2 



po\a 2 



The coefficient p reflects the correlation between the movements in the asset values of 
the two firms. Then the probability of firm i defaults at time t can be easily calculated 
as, 

X t (0)-ln(^ =9 N (_^ 

sTt, 



P(t) = 2-N 



where 



o l \ft 
Xi(0) - ln(«i) 



(3) 



is the standardized distance of firm i to its default point and N(-) denotes the cumula- 
tive probability distribution function for a standard normal variable. 

Furthermore, if we assume that jij = 7,;, then the probability of at least one firm 
defaults by time t can be written as dZhou[|2001al) : 



Piuj(t) 



1 



2r 



E- six 
n 



n=l,3,' 



rnrtto 
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-it 



Ii 



where I v (z) is the modified Bessel function I with order v and 

if p < 0, 
otherwise 



(4) 



tan 



p 



7r + tan 
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( ) 



7r + tan 1 ( Z z^/- P z ~j > otherwise, 
ro = Z 2 /sm(6 ). 
Then, the default correlation of these two firms is 

(t) = m) + p.jt) - p^p.jt) - p iU3 {t) 
^p t m-mmm-m)} ' 

However, none of the above known models includes jumps in the processes. At the 
same time, it is well-known that jumps are a major factor in the credit risk analysis. 
With jumps included in such analysis, a firm can default instantaneousl y because of a 



sudden drop in its value which is impossible under a diffusion process. iZhoul (1200 lbl) 
has shown the importance of jump risk in credit risk analysis of an obligor. He im- 
plemented a simulation method to show the effect of jump risk in the credit spread of 
defaultable bonds. He showed that the misspecification of stochastic processes gov- 
erning the dynamics of firm value, i.e., falsely specifying a jump-diffusion process as 
a continuous Brownian motion process, can substantially understate the credit spreads 
of corporate bonds. Therefore, for multiple processes, considering the simultaneous 
jumps can be a better way to estimate the correlated default rates. The multivariate 
jump-diffusion processes can provide a convenient way to describe multivariate and 
correlated processes with jumps. 

2.2 Multivariate jump-diffusion processes 

Let us consider a complete probability space (fi, F, P) with information filtration (Ft). 
Suppose that X t = ln(Vt) is a Ma rkov process in som e state space D C W 1 , solving 
the stochastic differential equation dDuffie etalLbOQOl) : 

dX t = fi(X t )dt + a(X t )dW t + dZ u (6) 

where W is an (F t )-standard Brownian motion in K"; /j, : D — > W\ a : D — > R raxn , 
and Z is a pure jump process whose jumps have a fixed probability distribution v on 
l n such that they arrive with intensity {\(X t ) : t > 0}, for some A : D — > [0, oo 



Under these cond i tions, the above model is reduced to an affine model if (lAhangarani , 
2005HDuffieeTall2000h : 



li(X t ,t) = K +K 1 Xt 

(a(X u t)a(X u t) T ) lJ = (Hoh + (H^X, 

\(X t ) = l +h-X t , (7) 

where K = (Kq^Kx) e M" x R nxn , H = (i?o,-Hi) G K" x ™ X R nxnxn , I = 
(l ,h) G R n x R nxn . 

As we mentioned, one of the major problems in the credit risk analysis is to estimate 
the default rate of a firm during a given time horizon. This problem is reduced to a first 
passage time problem. In order to obtain a computable multi-dimensional formulas of 
FPT distribution, we need to simplify Eq. (O and (|7]i as follows, 
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1 . Each Wt in Eq. (|6]l is independent; 



2. K\ = 0, Hi = and Zi = in Eq. (0 which means the drift term, the diffusion 
process (Brownian motion) and the arrival intensity are independent of the state 
vector X t ; 

3. The distribution of jump-size Z t is also independent with respect to X t . 
In this scenario, we can rewrite Eq. (|6]l as 

dX t = fidt + adW t + dZ t , (8) 

where 

fx = Kq, cra T = H , A = l . 



2.3 First passage time distribution of Brownian bridge 

Although for jump-diffusion processes, the closed form solutions are usually unavail- 
able, yet between each two jumps the process is a Brownian bridge for univariate jump- 
diffusion process. lAtiva et a have deduced one-dimensional first passage time 



distribution in time horizon [0, T]. In order to evaluate multiple processes, we obtain 
multi-dimensional formulas from Eq. <[8j and reduce them to computable forms. 

First, let us consider a firm i, as described by Eq. dS), such that its state vector Xi 
satisfies the following SDE: 

dXi = Hidt + UijdWj + dZi 
j 

= mdt + aidWi + dZi, (9) 
where Wi is a standard Brownian motion and erj is: 



We assume that in the interval [0, T], total number of jumps for firm i is Mj. Let the 
jump instants be Ti, T^, • • • , 7W 4 . Let Tq = and Tmi+i = T. The quantities r ? equa l 
to interjump times, which is Tj — Tj-i. Following the notation of lAtiva et a 3 d2005l) . 



let Xi(Tj~) be the process value immediately before the jth jump, and Xi{T^~) be the 
process value immediately after the jth jump. The jump-size is Xi(Tj) — Xi(T~), 
and we can use such jump-sizes to generate Xj(T^) sequentially. 

If we define Ai(t) as the event consisting of process crossed the threshold level 
Di(t) for the first time in the interval [t, t + dt], then we have 

g ij (t)=p(A i (t) g dtlXifT+^XiiTf)). (10) 

If we only consider one interval [T}-i , Tj], we can obtain 

XAT? ,) - DAt) 3 

9«W - J % , ^ (t-Tj^-iiTj-tyi 
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* cxp ( WF^I J 

* cxp ( n-T^of J' 

where 

[^(r+ 1 )-x t (r J -) + ^r J ] 2 ^ 



: exp 



2 

3 w i 



After getting result in one interval, we combine the results to obtain the density for 
the whole interval [0,T]. Let B(s) be a Brownian bridge in the interval [Tj_ x ,Tj] with 
B{Tf_ x ) = Xi(Tf_ x ) and B{Tr) = X t (Tr). Then the probability that the minimum 
of B(si) is always above the boundary level is 



Pij = L ^^B(a i )>D i (t)\B(T+ 1 )=X i (T+ 1 ),B(Tf)=X i (Tf 

\ 3 — 




2[Xj (T+_ t ) - J>, (t )] [X, (T- ) - D z (t)] 
7^3 



if X i (Tf)> A(i),, 12) 
otherwise. 



This implies that £?(s,) is below the threshold level, which means the default hap- 
pens or already happened, and its probability is 1 — P^. Let L(si) = Li denote 
the index of the interjump period in which the time Si (first passage time) falls in 
p£ 4 _i, IlJ. Also, let represent the index of the first jump, which happened in the 
simulated jump instant, 

h = min(j : X, (T fc ~) > A (t); fc = 1, . . . , j, and 

Xi{T+) > Di(t);k= 1, ... ,j - 1, and ^(T+) < £>i(t)).(13) 

If no such Ii exists, then we set /; = 0. 

By combining Eq. (fTTT i. (fl"2l) and (1131) . we get the probability of crossing the 
boundary level in the whole interval [0, T] as 

P(M*i) e rfslXi^t ^XiiTpJ = l,...,Mi + 1) 

ffji, (si) llfcii 1 ^»fe if -^t < ^ or ^ = 0) 

ftii(si) n^ 1 + uiu p l kS{ Sl - r h ) if u = i h d4) 

if L % > h, 

where S is the Dirac's delta function. 



2.4 The kernel estimator 

For firm i, after generating a series of first passage times s;, we use a kernel density 
estimator with Gaussian kernel to estimate the first passage time density (FPTD) /. 
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The kernel density estimator is based on centering a kernel function of a bandwidth as 
follows: 

f=-=Y l K(h,t-8i), (15) 



N 

i=l 



where 

K(h, t— Si) = — -^=- exp ' 



T he optimal bandwidth in the kernel function K can be calculated as ( Silverman! 
1983) : 



-0.2 

F"\2, 



= 2iVy^ / (/ t ") 2 d* , (16) 



where N is the number of generated points and f t is the true density. H ere we use 
the app roximation for the distribution as a gamma distribution, proposed bv lAtiva et al 
d2Q05h : 

ft = ^- ) t fj - 1 cxp(-at). (17) 



So the integral in Eq. ( fT6l > becomes: 



/ (/n 2 dt = j2 



W iai r(2/3 - i) 



(18) 



^2^-) ( r(/3))2' 

where 

Wi = A 2 , W 2 =2AB, W 3 = B 2 + 2AC, Wi = 2BC, W 5 = C 2 , 

and 

A = a 2 , B = -2a(p-l), C = (0 - l)(/3 - 2). 

From Eq. dT8b . it follows that in order to get a nonzero bandwidth, we have to have 
constraint (3 to be at least equal to 3. 

The kernel estimator for the multivariate case involves the evaluation of joint condi- 
tional interjump first passage time density, as discussed in Section[3] The methodology 
for such an evaluation is quite involved compared to the one-dimensional case and we 
will focus on these details elsewhere. In what follows we highlight the main steps of 
the procedure. 



3 The methodology of solution 

First, let us recall the conventional Monte-Carlo procedure in application to the anal- 
ysis of the evolution of firm Xi within the time horizon [0, T]. We divide the time 
horizon into n small intervals [0, t{\, [ti, £2], ■ • •, [t n -i,T] as displayed in Fig. 01 a )- 
In each Monte Carlo run, we need to calculate the value of Xi at each discretized 
time t. As usual, in order to exclude discretization bias, the number n must be large. 
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This procedure exhibits substantial computational difficulties when applied to jump- 
diffusion processes. Indeed, for a typical jump-diffusion process, as shown in Fig. 
[TJa), let Tj_i and Tj be any successive jump instants, as described above. Then, in the 
conventional Monte Carlo method, although there is no jump occurring in the interval 
[Tj-i, Tj], yet we need to evaluate at each discretized time t in [T 3 _i, Tj]. This 
very time-consuming procedure results in a serious shortcoming of the conventional 
Monte-Carlo methodology. 

Conventional Monte-Carlo method 



(a) 

Uniform sampling method 
Sj (uniform distribution) 



X f T + \ X ( T~ 1 

( \ ' 'K Brownian bridge $H i[ - ' ' 

(b) 

Figure 1 : Schematic diagram of (a) conventional Monte Carlo and (b) uniform sam- 
pling (UNIF) method. 

To remedy the situation, two modifications of the con ventional procedure were 



recently proposed dAtiva et all I2005t iMetwallv et all |2002) that allow us a potential 



speed-up of the conventional methodology in 10-30 times. One of the modifications, 
the uniform sampling method, involves samplings using the uniform distribution. The 
other is the inverse Gaussian density sampling method. Both methodologies were de- 
veloped for the univariate case. 

The major improvement of the uniform sampling method is based on the fact that 
it only evaluates X, at generated jump times, while between each two jumps the pro- 
cess is a Brownian bridge (see Fig. |TJb)). Hence, we just consider the probability of 
Xi crossing the threshold in (Tj-i,Tj) instead of evaluating X; at each discretized 
time t. More precisely, in the uniform sampling method, we assume that the values 
of Xj(Tj 1 l 1 ) and Xj(T~) are known as two end points of the Brownian bridge, the 
probability of firm i defaults in (Tj_i,Tj) is 1— Py which can be computed according 
to Eq. (TT2b . Then we generate a variable from a distribution uniform in the inter- 
val [Tj-i, Tj-i H — \~lp.~ 1 ]■ If the generated point Si falls in the interjump interval 
[Tj_i, Tj], then we have successfully generated the first passage time s; and can ne- 
glect the other intervals and perform another Monte Carlo run. On the other hand, if the 
generated point Si falls outside the interval [Tj_i , Tj] (which happens with probability 
Pij), then that point is "rejected". This means that no boundary crossing has occurred 
in the interval, and we proceed to the next interval and repeat the whole process again. 

Note that the generated Si is not obtained according to conditional boundary cross- 
ing densi ty gn(sj) a s descr ibed by Eq. ( fTTT i. In order to obtain an appropriate density 



estimate, lAtiva et all ( 120051) proposed that the right hand side summation in Eq. (T5[ 
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can be viewed as a finite sample estimate of the following: 

E gii(si)[ K {h,t- Si)] = / g tj {si)K(h,t- Si)dsi 



T 3 -i 

Tj T±-x 
1 - Pio 



Eu( Si )[9ij(si)K(h,t- Si )], (19) 



where E g ..^ means the expectation of Sj, where Sj obeys the density gij(si). U (s.;) 
is the uniform density in 2j_i + -j^-fr-^-] from which we sample the point s;. 

Therefore, we should weight the kernel with ( "^Zp.T 1 ) 9ij( s i) to obtain an estimate 
for the true density. 

For the multidimensional density estimate we need to evaluate the joint conditional 
boundary crossing density. This problem can be divided into seve ral one-dimensiona l 
density estimate subproblems if the processes are non-correlated dZhang et al , 120061) . 



As for the multivariate correlated processes, the joint density becomes very compli- 
cated and there are usually no ana lytical solutions for higher-dimensional processes 



dSong et all l2006t IWise et al 120041) . We will not consider this problem in the current 
contribution. 

Instead, we focus on the further development of the uniform sampling (UNIF) 
method and extend it to multivariate and correlated jump-diffusion processes. In or- 
der to implement the UNIF method for our multivariate model as described in Eq. ©, 
we need to consider several points: 

1 . We assume that the arrival rate A for the Poisson jump process and the distribu- 
tion of (Tj —Tj—i) are the same for each firm. As for the jump-size, we generate 
them by a given distribution which can be different for different firms to reflect 
specifics of the jump process for each firm. 

2. We exemplify our description by considering an exponential distribution (mean 
value Ht) for (Tj — Tj_i) and a normal distribution (mean value \ij and stan- 
dard deviation aj) for the jump-size. We can use any other distribution when 
appropriate. 

3. An array IsDefault (whose size is the number of firms denoted by Afg rm ) is 
used to indicate whether firm % has defaulted in this Monte Carlo run. If the firm 
defaults, then we set IsDefault(i) = 1, and will not evaluate it during this 
Monte Carlo run. 

4. Most importantly, as we have mentioned before, the default events of firm i 
are inevitably correlated with other firms, for example firm i + 1. The default 
correlation of firms i and i + 1 is described by Eq. @- Hence, firm i's first 
passage time Si is indeed correlated with s^+i - the first passage time of firm 
i + 1. We must generate several correlated s; in each interval [Tj_i, + 
^"p -1 ] which is the key point for multivariate correlated processes. 

Note that the assumption based on using the same arrival rate A and distribution of 
(Tj — Tj-i) for different firms may seem to be quite idealized. One may argue that the 
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arrival rate A for the Poisson jump process should be different for different firms, which 
implies that different firms endure different jump rates. However, if we consider the 
real market economy, once a firm (called firm "A") encounter sudden economic hazard, 
its correlated firms may also endure the same hazard. Furthermore, it is common that 
other firms will help firm "A" to pull out, which may result in a simultaneous jump 
for them. Therefore, as a first step, it is reasonable to employ the simultaneous jumps' 
processes for all the different firms. 

Next, we will give a brief description of the sum-of-uniforms method which is used 
to generate correlated uniform random variables, followed by the description of the 
multivariate and correlated UNIF method and the model calibration. 



3.1 Sum-of-uniforms method 



In the above sections, we have reduced the solution of the original problem to a series 
of one-dimensional jump-diffusion processes as described by Eq. (|9]). The first passage 
time distribution in an interval [Tj_i, Tj] (between two successive jumps) was obtained 
in section |2~3l As mentioned, the default events of firm i are inevitably correlated with 
other firms, for example firm In this contribution, we approximate the correlation 
of Si and Si+i as the default correlation of firm si and i + 1 by the following formula: 



p(si, Sj + i) w p hl+1 (t) 



Pi(t) + p i+1 (t) - Pj(t)P i+1 (t) - iWi(*) 



(20) 



where t can be chosen as the midpoint of this interval. 

Therefore, we need to generate several correlated Sj in [Tj_i,Tj_x 
whose correlations can be described by Eq. d20l i. Let us introduce a new variable 



Tj-Tj. 

l-Pi: 



T- — T- 
± j ± j- 

1-Pi. 



-, then we have s, = bijYi +Tj—i, where Yi are uniformly distributed in 



[0, 1]. Moreover, the correlation of and Yi + i is given by p(s;, Sj+i). 

Now we can generate the correla ted uniform random variables Y\ , Y % , ■ ■ ■ by using 
the sum-of-uniforms (SOU) method dChenj|2005llWillemain et all 1 1993b in the follow- 
ing steps: 

1. Generate Y\ from numbers uniformly distributed in [0, 1]. 

2. For i = 2,3, • • generate Wi ~ U(0, d.--\ a), wher e U (0. Cj_i j) denotes a 
uniform random number over range (0,Ci_i,i). IChenl J2005h has obtained the 
relationship of parameter Cj_i j and the correlation p(si-i, s<) (abbreviated as 
Pi-i,i) as follows: 



Correlation Ci-n > 1 



-^-{Pi-i,i > -o.7) 



Ci_l,i < 1 



1 - 0.5^ + O^Lite-i,* > 0.7) 
-1 + 0.5c?_! 4 - Q.24_ x A Pi -iA < -0.7) 



Pi-1A > 
Pi-IA < 



If Yi-i and Y{ are positively correlated, then let 

Z i =Y i - 1 + W i . 
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If Yi-i and Y t are negatively correlated, then let 



Z i = l-Y i _ 1 +W i . 
Let Yi = F(Zi), where for Cj_i,i > 1, 

f Z 2 /(2c 8 _ M ), 0<Z<1, 
F(Z)=^ (2Z - l)/(2c i _i >i ), l^^^Ci-i,,, 

[ l-(l + c t _i, t -Z)7(2 Cl „ M ), Ci_ M < Z < 1 + Ci_ M , 

and for < Cj_x,j < 1, 

f Z 2 /(2c 4 _ M ), 0<Z<Ci_i,i, 
F(Z) = <^ (2Z - c i _ 1> 0/2, Q_ M <Z<1, 

[ l-(l + c l _ 1 . l -Z) 2 /(2c l _ 1 , l ), l<Z<l + q_ M . 

3.2 Uniform sampling method 

In this subsection, we will describe our algorithm for multivariate jump-diffusion pro- 
cesses, whic h is an extension o f the one-dimensiona l case developed earlier by other 
authors (e.g. lAtiva etail d2005l) : lMetwallv et all d2002l) ). 

Consider N& xm firms in the given time horizon [0, T]. First, we generate the jump 
instantly by generating interjump times (Tj —Tj_i) and set all the IsDef ault(i) = 
Q(i = 1, 2, • • • , Nfi rm ) to indicate that no firm defaults at first. 

From Fig. [TJb) and Eq. (0, we can conclude that for each process Xi we can make 
the following observations: 

1. If no jump occurs, as described by Eq. (0, the interjump size (Xi(T^) — 
^(X^j)) follows a normal distribution of mean Hi(Tj — T}_i) and standard 
deviation Oi \jTj — Tj~\_. We get 



X t {T7) ~ XiiTf^) + in (Tj - T^) + aiN(0, T 3 - T 3 ^) 

iV firm 



k=l 

where the initial state is Xi(0) = Xj(Tg"). 

2. If jump occurs, we simulate the jump-size by a normal distribution or another 
distribution when appropriate, and compute the postjump value: 

Xi(T+) = XiiTf) + Z^). 

This completes the procedure for generating beforejump and postjump values Xi (T~ ) 
and Xi(Tj). As before, j = 1, • • • , M where M is the total number of jumps for all 
the firms. We compute Pjj according to Eq. ( fT2] >. To recur the first passage time den- 
sity (FPTD) fi(t), we have to consider three possible cases that may occur for each 
non-default firm i: 
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1. First passage happens inside the interval. We know that if X^T^) > 
Di(Tj^i) and Xi(Tj~) < Di(Tj), then the first passage happened in the time 
interval [Tj-±, Tj]. To evaluate when the first passage happened, we introduce 
a new viable by as fry = "^Jp.T 1 ■ We generate several correlated uniform 
numbers Yi by using the SOU method as described in Section [3~Tl then compute 
Si = bijYi + Tj-i. If Si belongs to interval [Tj_i,Tj], then the first passage 
time occurred in this interval. We set IsDef ault(i) = 1 to indicate firm i 
has defaulted and compute the conditional boundary crossing density gij(si) ac- 
cording to Eq. ( fTTT i. To get the density for the entire interval [0, T], we use 

fi,n(t) = ( ^p.T 1 ) 9ij( s i) * K(h optl t — Si), where n is the iteration number 
of the Monte Carlo cycle. 

2. First passage does not happen in this interval. If doesn't belong to interval 
[Tj-i, Tj], then the first passage time has not yet occurred in this interval. 

3. First passage happens at the right boundary of the interval. If Xi(Tf) < 

Di (Tj) and Xj(Tr) > D^) (see Eq. (TBI)'), then T u is the first passage time 

and Li = j, we evaluate the density function using kernel function fi. n {t) = 
K(hopt,t- TjJ, and set IsDefault(i) = 1. 

Next, we increase j and examine the next interval and analyze the above three cases 
for each non-default firm again. After running N times Monte Carlo cycle, we get the 
FPTD of firm i as /;(*) = % £f =1 f hn {t). 



3.3 Model calibration 

We need to calibrate the developed model, in other words, to numerically choose or 
optimize the parameters, such as drift, volatility and jumps to fit the most liquid market 
data. This can be done by applying the least-square method, minimizing the root mean 
square error {rinse) given by: 



(Market price — Model price) 5 




Number of derivatives 



Luciano et all (2005) have used a set of European call options C(k, T) as their model 



price to calibrate their model parameters. 

However, as demonstrated in Section |4] for a number of practically interesting 
cases, there is no option value that can be used to calibrate our model, so we have 
to use the historical default data to optimize the parameters in the model. As men- 
tioned in Sections 12.41 and 13.21 after Monte Carlo simulation we obtain the estimated 
density fi(t) by using the kernel estimator method. The cumulative default rates for 
firm i in our model is defined as, 

Pi(t) = I fi(r)dr. (21) 
Jo 
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Then we minimize the difference between our model and historical default data 
Ai (t) to obtain the optimized parameters in the model (such as er^ , arrival intensity A 
inEq. ©): 

/ i ; 



4 Applications and discussion 

In this section, we demonstrate the developed model at work for analyzing the default 
events of multiple correlated firms via a set of historical default data. 

4.1 Density function and default rate 

First, for completeness, le t us consider a set of historical default data of differently 



rated firms as presented by lZhoul (1200 lal) . Our first task is to describe the first passage 
time density functions and default rates of these firms. 

Since there is no option value that can be used, we will employ Eq. (l22l) to optimize 
the parameters in our model. For convenience, we reduce the number of optimizing 
parameters by: 

1. Setting X(0) = 2 and = 0. 

2. Setting the g rowth rate 7 of debt value equivalent to the growth rate /1 of the 



firm's value dZhoul 1200 lal) . so the default of firm is non-sensitive to /i. In our 



computations, we set /j = —0.001. 

3. The interjump times (Tj — Tj-i) satisfy an exponential distribution with mean 
value equals to 1 . 

4. The arrival rate for jumps satisfies the Poisson distribution with intensity param- 
eter A, where the jump size is a normal distribution Z t ~ N(fj,z, <Jz)- 

As a result, we only need to optimize tr, A, ^z, uz for each firm. This is done 
by minimizing the differences between our simulated default rates and historical data. 
Moreover, as mentioned above, we will use the same arrival rate A and distribution of 
(Tj — Tj_i) for differently rated firms, so we first optimize four parameters for, e.g., 
the A-rated firm, and then set the parameter A of other three firms the same as A's. 

The minimization was performed by the using quasi-Newton procedure imple- 
mented as a Scilab program. The optimized parameters for each firm are described 
in Table □ 

By using these optimized parameters, we carried out the final simulation with 
Monte Carlo runs N = 500, 000. The estimated first passage time density function 
of these four firms are shown in Fig. |2] The simulated cumulative default rates (line) 
together with historical data (squares) are given in Fig. [3] The theoretical data d enoted 



as circ les in Fig. [3]were computed by using Eq. (01 where Zi were evaluated in ( iZhou , 



2001 a) as 8.06, 6.46, 3.73 and 2.10 for A-, Baa-, Ba- and B-rated firms, respectively. 
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Table 1 : Optimized parameters for differently rated firms by using the UNIF method. 
The optimization was performed by using the quasi-Newton procedure implemented 
as a Scilab program. In each step of the optimization, we choose the Monte Carlo runs 
iV = 50, 000. 





a 


A 




&z 


A 


0.0900 


0.1000 


-0.2000 


0.5000 


Baa 


0.0894 


0.1000 


-0.2960 


0.6039 


Ba 


0.1587 


0.1000 


-0.5515 


1.6412 


B 


0.4500 


0.1000 


-0.8000 


1.5000 



In Table [2] we give the optimal bandwidth and parameters a, (3 for the true density 
estimate. 




1 3 5 7 9 11 13 15 17 19 1 3 5 7 9 1 1 13 15 17 19 

Year Year 



Figure 2: Estimated density function for differently rated firms. All the simulations 
were performed with Monte Carlo runs N = 500, 000. 

Based on these results, we conclude that: 

1. Simulations give similar or better results to the analytical results predicted by 
Eq. ©. 

2. A- and Baa-rated firms have a smaller Brownian motion part. Their parameters 
a are much smaller than those of Ba- and B-rated firms. 

3. The optimized parameters a of A- and Baa-rated firms are similar, but the jump 
parts (fiz,&z) are different, which explains their different cumulative default 
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Figure 3: Historical (squares), theoretical (circles) and simulated (line) cumulative 
default rates for differently rated firms. All the simulations were performed with Monte 
Carlo runs N = 500,000. 



Table 2: The optimal bandwidth h op t, parameters a, f3 for the true density estimate 
of differently rated firms. All the simulations were performed with Monte Carlo runs 
N = 500, 000. 





a 


P 


Optimal bandwidth 


A 


0.206699 


3 


0.655522 


Baa 


0.219790 


3 


0.537277 


Ba 


0.252318 


3 


0.382729 


B 


0.327753 


3 


0.264402 
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rates and density functions. Indeed, Baa-rated firm may encounter more severe 
economic hazard (large jump-size) than A-rated firm. 

4. As for Ba- and B-rated firms, except for the large a, both of them have large 
/iz and especially large az, which indicate that the loss due to sudden economic 
hazard may fluctuate a lot for these firms. Hence, the large a, fj,z and az account 
for their high default rates and low credit qualities. 

5. From Fig. |2j we can conclude that the density functions of A- and Baa-rated 
firms still have the trend to increase, which means the default rates of A- and 
Baa-rated firms may increase little faster in future. As for Ba- and B-rated firms, 
their density functions have decreased, so their default rates may increase very 
slowly or be kept at a constant level. Mathematically speaking, the cumulative 
default rates of A- and Baa-rated firms are convex function, while the cumulative 
default rates of Ba- and B-rated firms are concave. 

4.2 Correlated default 

Our final example concerns with the default correlation of two firms. If we do not 
include jumps in the model, the default correlation can be calculated by using Eq. (f3), 
(O and ©. In Tables 06] we pre sent compariso ns of our results with those based on 
closed form solutions provided bv lZhoul (1200 lal) with p = 0.4. 

Table 3: One year default correlations (%). All the simulations are performed with the 

Monte Carlo runs N = 500,000 

UNIF Zhou (2001a) 





A 


Baa 


Ba 


B 


A 


Baa 


Ba 


B 


A 


-0.01 








0.00 








Baa 


-0.02 


3.69 






0.00 


0.00 






Ba 


2.37 


4.95 


19.75 




0.00 


0.01 


1.32 




B 


2.80 


6.63 


22.57 


26.40 


0.00 


0.00 


2.47 


12.46 



Table 4: Two year default correlations (%). All the simulations are performed with the 

Monte Carlo runs N = 500,000 

UNIF Zhou (2001a) 





A 


Baa 


Ba 


B 


A 


Baa 


Ba 


B 


A 


2.35 








0.02 








Baa 


2.32 


4.25 






0.05 


0.25 






Ba 


4.17 


7.17 


20.28 




0.05 


0.63 


6.96 




B 


4.73 


8.23 


23.99 


29.00 


0.02 


0.41 


9.24 


19.61 



Next, let us consider the default correlations under the multivariate jump-diffusion 
processes. We use the following conditions in our multivariate UNIF method: 

1. Setting X(0) = 2 and ln( K ) = for all firms. 
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Table 5: Five year default correlations (%). All the simulations are performed with the 

Monte Carlo runs N = 500,000 

UNIF Zhou (2001a) 





A 


Baa 


Ba 


B 


A 


Baa 


Ba 


B 


A 


6.45 








1.65 








Baa 


6.71 


9.24 






2.60 


5.01 






Ba 


7.29 


10.88 


22.91 




2.74 


7.20 


17.56 




B 


6.77 


10.93 


22.97 


27.93 


1.88 


5.67 


18.43 


24.01 



Table 6: Ten year default correlations (%). All the simulations are performed with the 
Monte Carlo runs N = 500, 000 



UNIF Zhou (2001a) 





A 


Baa 


Ba 


B 


A 


Baa 


Ba 


B 


A 


8.79 








7.75 








Baa 


10.51 


13.80 






9.63 


13.12 






Ba 


9.87 


14.23 


22.50 




9.48 


14.98 


22.51 




B 


8.50 


12.54 


20.49 


24.98 


7.21 


12.28 


21.80 


24.37 



2. Setting 7 = \i and /i = —0.001 for all firms. 

3. Since we are considering two correlated firms, we choose cr as, 



Oil 
021 



012 
022 



(23) 



where aa T = Hn such that, 



00 1 =H = 



P12C1C2 



P12CT1C2 



and 



a 



'11 
2 

21 



a 



'12' 
2 

22- 



Pl2 



011021 + 012022 
0102 



(24) 



In Eq. d24l i. p\2 reflects the correlation of diffusion parts of the state vectors of 
the two firms. In order to compare with the standard Brownian motion and to 
eval uate the defaul t correlations between different firms, we set all the pi2 = 0.4 
as in lZhoul(l2001al) . Furthermore, we use the optimized o\ and 02 in TableQ]for 
firm 1 and 2, respectively. Assuming a\i = 0, we get, 

011 =01, 

012 = 0, 
021 = P1202, 



022 = yl 



P? 2 02- 
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4. The arrival rate for jumps satisfies the Poisson distribution with intensity pa- 
rameter A = 0.1 for all firms. The jump size is a normal distribution Z t ~ 
N(/iZi, fz,), where and az t can be different for different firms to reflect 
specifics of the jump process for each firm. We adopt the optimized parameters 
given in Table Q] 

5. As before, we generate the same interjump times (Tj — Tj-i) that satisfy an 
exponential distribution with mean value equals to 1 for each two firms. 

We carry out the UNIF method to evaluate the default correlations via the following 
formula: 

N t^i V P hn(t)(l - P ltn (t))P 2 , n (t)(l - P 2>n (t)) ' 

where P\2, n {t) is the probability of joint default for firms 1 and 2 in each Monte Carlo 
cycle, Pi^ n (t) and P^.n (t) are the cumulative default rates of firm 1 and 2, respectively, 
in each Monte Carlo cycle. 

The simulated default correlations for one-, two-, five- and ten-years are given 
in Table 06] All the simulations were performed with the Monte Carlo runs N = 
500, 000. Comparing those simulated default correlations with the theoretical data for 
standard Brownian motions, we can conclude that 



1 . Similarly to conclusions of lZhoul (1200 lal) , the default correlations of same rated 



firms are usually large compared to differently rated firms. Furthermore, the 
default correlations tend to increase over long horizons and may converge to a 
stable value. 

In our simulations, the one year default correlations of (A,A) and (A,Baa) are 
negative. This is because they seldom default jointly during one year. Note, how- 
ever, that the def ault correlatio ns of other firms are positive and usually larger 



than the results in lZhoul (12001 ah . 



3. For two and five years, the default correlations of different firms increase. This 
can be explained by the fact that their individual first passage time density func- 
tions increase during these time horizon, hence the probability of joint default 
increases. 

4. As for ten year default correlations, our simulated results are almost identical to 
the theoretical data for standard Brownian motions. The differences are that the 
default correlations of (Ba,Ba), (Ba,B) and (B,B) decrease from the fifth year to 
tenth year in our simulations. The reason is that the first passage time density 
function of Ba- and B-rated firms begin to decrease from the fifth year, hence the 
probability of joint default may increase slowly. 

5 Conclusion 

In this contribution, we have analyzed the credit risk problems of multiple correlated 
firms in the structural model framework, where we incorporated jumps to reflect the ex- 
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ternal shocks or other unpredicted events. By combining the fast Monte-Carlo method 
for one-dimensional jump-diffusion processes and the generation of correlated multidi- 
mensional variates, we have developed a fast Monte-Carlo type procedure for the anal- 
ysis of multivariate and correlated jump-diffusion processes. The developed approach 
generalizes previously discussed non-correlated jump-diffusion cases for multivariate 
and correlated jump-diffusion processes. Finally, we have applied the developed tech- 
nique to analyze the default events of multiple correlated firms via a set of historical 
default data. The developed methodology provides an efficient computational tech- 
nique that is applicable in other areas of credit risk and pricing options. 
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